body .main-container {
max-width: 100% !important;
width: 100% !important;
}
body {
max-width: 100% !important;
}
body, td {
font-size: 16px;
}
code.r{
font-size: 14px;
}
pre {
font-size: 14px
}knitr::opts_chunk$set(message = FALSE,
warning = FALSE,
echo = TRUE,
include = TRUE,
fig.align = "center",
out.width = "100%"
)
set.seed(1982)library(INLA)
library(here)data(Seeds)
df <- data.frame(y = Seeds$r, Ntrials = Seeds$n, Seeds[, 3:5])
df |> head(2) |> knitr::kable(caption = "Data")| y | Ntrials | x1 | x2 | plate |
|---|---|---|---|---|
| 10 | 39 | 0 | 0 | 1 |
| 23 | 62 | 0 | 0 | 2 |
The model specification is
\[ y_i\sim\text{Binomial}\left(N_i, p_i\right)\label{eq1}\tag{A} \]
\[ p_i = \text{invlogit}(\eta_i) = \frac{\exp(\eta_i)}{1+\exp(\eta_i)} = \frac{1}{1+\exp(-\eta_i)}\qquad \left(\text{i.e., }\eta_i = \text{logit}(p_i) =\log\left(\frac{p_i}{1-p_i}\right) \right) \]
\[ \eta_i = \beta_0 + \beta_1x_{1i} + \beta_2x_{2i} +v_i \]
\[ \beta_i \sim \text{N}(0, 1000) \]
\[ v_i\sim\text{N}(0,\sigma^2) \] \[ \sigma\sim\text{Exp}(\lambda) \qquad \left(\text{i.e., }\pi(\sigma) = \lambda e^{-\lambda\sigma}\right) \]
\[ \texttt{inla()}\text{ deals with }\sigma^{-2} = \tau \sim\pi(\tau) =\pi(\sigma) \left|\dfrac{d\sigma}{d\tau}\right|= \dfrac{\lambda}{2}\tau^{-3/2}e^{-\lambda \tau^{-1/2}} \]
\[ \text{internally }\texttt{inla()}\text{ deals with }\log(\tau) = \theta \sim\pi(\theta) = \pi(\tau)\left|\dfrac{d\tau}{d\theta}\right| = \dfrac{\lambda}{2}\exp\left(-\dfrac{\theta}{2} - \lambda \exp\left(-\dfrac{\theta}{2}\right)\right) \]
\[ \lambda \text{ is such that } \pi(\sigma>1) = 0.01\qquad\left(\text{i.e., }0.01 = \pi(\sigma>1) = 1- \pi(\sigma\leq1) = e^{-\lambda}\implies \lambda = -\dfrac{\log(0.01)}{1}\approx4.6\right)\label{eq2}\tag{B} \]
# call to inla
hyper <- list(theta = list(prior = "pc.prec", param = c(1,0.01)))
formula1 = y ~ x1 + x2 + f(plate,
model = "iid",
hyper = hyper)
model1 = inla(formula = formula1,
data = df,
family = "binomial",
Ntrials = Ntrials,
control.family = list(control.link = list(model = "logit")),
control.predictor = list(compute = TRUE, link = 1))
summary(model1)## Time used:
## Pre = 0.373, Running = 0.159, Post = 0.0129, Total = 0.545
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## (Intercept) -0.387 0.181 -0.739 -0.390 -0.017 -0.390 0
## x1 -0.360 0.231 -0.838 -0.353 0.076 -0.353 0
## x2 1.033 0.221 0.589 1.034 1.469 1.034 0
##
## Random effects:
## Name Model
## plate IID model
##
## Model hyperparameters:
## mean sd 0.025quant 0.5quant 0.975quant mode
## Precision for plate 19.78 32.61 2.96 10.55 84.94 6.06
##
## Marginal log-Likelihood: -68.44
## is computed
## Posterior summaries for the linear predictor and the fitted values are computed
## (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')# data generation
N = 3600
sig.epsilon = 0.5
sig.u = 1
sig.seasonal = 2
rho = 0.90
u = arima.sim(list(order = c(1, 0 ,0), ar = rho), n = N,sd = 1)
u = u/sd(u)*sig.u
seas.coeff = (0:11)*(1:12-12)
seas = rep(seas.coeff, N/12)
seas = drop(scale(seas))*sig.seasonal
gaussian.error = rnorm(N, 0, sd = sig.epsilon)(prec.sig.epsilon = sig.epsilon^(-2))## [1] 4(prec.sig.u = sig.u^(-2))## [1] 1(prec.sig.seasonal = sig.seasonal^(-2))## [1] 0.25# the model
y = u + seas + gaussian.errordf = data.frame(y = y, t = 1:N, year = rep(1:12, N/12))
df |> head(2) |> knitr::kable(caption = "Data")| y | t | year |
|---|---|---|
| 1.9358080 | 1 | 1 |
| 0.8643198 | 2 | 2 |
\[ y_i = \beta_0 + u(t_i) + s_i + \epsilon_i,\; i =1,\dots, n \]
\[ u(t_1) \sim \text{N}(0, (\tau_u(1-\rho^2))^{-1}), \qquad u(t_i) = \rho u(t_{i-1}) + e_i, \qquad e_i \sim \text{N}(0, \tau_u^{-1} = \sigma_u^2) \qquad i =2,\dots, n, \qquad |\rho|<1\qquad \text{or}\qquad u(t_i) - \rho u(t_{i-1}) \sim \text{N}(0, \tau_u^{-1} = \sigma_u^2) \]
\[ \beta_0 \sim \text{N}(0, 1000) \]
\[ \epsilon_i \sim \text{N}(0, \tau_\epsilon^{-1} = \sigma_\epsilon^2) \]
\[ s_i-2s_{i+1}+s_{i+2}\sim\text{N}(0, \tau_s^{-1} = \sigma_s^2) \]
\[ \sigma_u^{-2} = \tau_u \sim\pi(\tau_u) = \dfrac{\lambda}{2}\tau_u^{-3/2}e^{-\lambda \tau_u^{-1/2}},\qquad\text{where }\lambda = -\dfrac{\log(0.5)}{0.1}\text{ is such that }\pi(\sigma_u>0.1) = 0.5 \]
\[ \rho\sim\pi(\rho) = \dfrac{\lambda\exp(-\lambda\sqrt{1-\rho})}{1-\exp(-\sqrt2\lambda)}\dfrac{1}{2\sqrt{1-\rho}},\qquad\text{where }\lambda,\text{ the solution of }\dfrac{\exp(-\lambda\sqrt{1-0.9})}{1-\exp(-\sqrt2\lambda)} = 0.5,\text{ is such that }\pi(\rho>0.9) = 0.5 \]
\[ \sigma_s^{-2} = \tau_s \sim\pi(\tau_s) = \dfrac{\lambda}{2}\tau_s^{-3/2}e^{-\lambda \tau_s^{-1/2}},\qquad\text{where }\lambda = -\dfrac{\log(0.5)}{0.1}\text{ is such that }\pi(\sigma_s>0.1) = 0.5 \]
\[ \sigma_\epsilon^{-2} = \tau_\epsilon \sim\pi(\tau_\epsilon) = \dfrac{\lambda}{2}\tau_\epsilon^{-3/2}e^{-\lambda \tau_\epsilon^{-1/2}},\qquad\text{where }\lambda = -\dfrac{\log(0.5)}{3}\text{ is such that }\pi(\sigma_\epsilon>3) = 0.5 \]
hyper.ar1 <- list(theta1 = list(prior="pc.prec", param = c(0.1, 0.5)), #sigma_u
theta2 = list(prior="pc.cor1", param = c(0.9, 0.5))) #rho
hyper.rw2 <- list(theta1 = list(prior="pc.prec", param = c(1.99, 0.99))) #sigma_s
hyper.family <- list(theta = list(prior="pc.prec", param = c(3, 0.5))) #sigma_epsilon
formula <- y ~ f(t, model = 'ar1', hyper = hyper.ar1) +
f(year, model = "rw2", hyper = hyper.rw2, cyclic = T, constr = T)
model2 <- inla(formula = formula,
data = df,
family = "gaussian",
control.predictor = list(compute = TRUE),
control.family = list(hyper = hyper.family))
summary(model2)## Time used:
## Pre = 0.192, Running = 1.1, Post = 0.164, Total = 1.45
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## (Intercept) 0.043 0.076 -0.106 0.043 0.192 0.043 0
##
## Random effects:
## Name Model
## t AR1 model
## year RW2 model
##
## Model hyperparameters:
## mean sd 0.025quant 0.5quant
## Precision for the Gaussian observations 3.990 0.184 3.640 3.986
## Precision for t 1.018 0.078 0.870 1.016
## Rho for t 0.907 0.009 0.889 0.908
## Precision for year 1.210 0.546 0.439 1.111
## 0.975quant mode
## Precision for the Gaussian observations 4.365 3.977
## Precision for t 1.178 1.013
## Rho for t 0.924 0.908
## Precision for year 2.542 0.929
##
## Marginal log-Likelihood: -4048.59
## is computed
## Posterior summaries for the linear predictor and the fitted values are computed
## (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')formula <- y ~ f(t, model = 'ar1') +
f(year, model = "rw2", cyclic = T, constr = T)
model2 <- inla(formula = formula,
data = df,
family = "gaussian",
control.predictor = list(compute = TRUE))
summary(model2)## Time used:
## Pre = 0.215, Running = 1.49, Post = 0.274, Total = 1.98
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## (Intercept) 0.043 0.076 -0.107 0.043 0.194 0.043 0
##
## Random effects:
## Name Model
## t AR1 model
## year RW2 model
##
## Model hyperparameters:
## mean sd 0.025quant 0.5quant
## Precision for the Gaussian observations 3.997 0.185 3.642 3.994
## Precision for t 1.007 0.078 0.859 1.005
## Rho for t 0.908 0.009 0.890 0.908
## Precision for year 1.570 0.633 0.651 1.464
## 0.975quant mode
## Precision for the Gaussian observations 4.371 3.990
## Precision for t 1.167 1.003
## Rho for t 0.925 0.908
## Precision for year 3.099 1.270
##
## Marginal log-Likelihood: -4057.16
## is computed
## Posterior summaries for the linear predictor and the fitted values are computed
## (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')# temp = read.csv("data/harmonised-unemployment-rates-mo.csv")
n = nrow(temp) - 1
df = data.frame(y = temp[1:n, 2], t = 1:n, dates = temp[1:n, 1])
df |> head(2) |> knitr::kable(caption = "Data")| y | t | dates |
|---|---|---|
| 9.1 | 1 | 2000-03 |
| 8.8 | 2 | 2000-04 |
\[ y_i = \beta_0+u(t_i)+\epsilon_i \]
\[ u(t_1) \sim \text{N}(0, (\tau_u(1-\rho^2))^{-1}), \qquad u(t_i) = \rho u(t_{i-1}) + e_i, \qquad e_i \sim \text{N}(0, \tau_u^{-1} = \sigma_u^2) \qquad i =2,\dots, n, \qquad |\rho|<1\qquad\text{or}\qquad u(t_i) - \rho u(t_{i-1}) \sim \text{N}(0, \tau_u^{-1} = \sigma_u^2) \]
\[ \beta_0 \sim \text{N}(0, 1000) \]
\[ \epsilon_i \sim \text{N}(0, \tau_\epsilon^{-1} = \sigma_\epsilon^2) \]
\[ \sigma_u^{-2} = \tau_u \sim\pi(\tau_u) = \dfrac{\lambda}{2}\tau_u^{-3/2}e^{-\lambda \tau_u^{-1/2}},\qquad\text{where }\lambda = -\dfrac{\log(0.5)}{0.02}\text{ is such that }\pi(\sigma_u>0.02) = 0.5 \]
\[ \rho\sim\pi(\rho) = \dfrac{\lambda\exp(-\lambda\sqrt{1-\rho})}{1-\exp(-\sqrt2\lambda)}\dfrac{1}{2\sqrt{1-\rho}},\qquad\text{where }\lambda,\text{ the solution of }\dfrac{\exp(-\lambda\sqrt{1-0.9})}{1-\exp(-\sqrt2\lambda)} = 0.5,\text{ is such that }\pi(\rho>0.9) = 0.5 \]
\[ \sigma_\epsilon^{-2} = \tau_\epsilon \sim\pi(\tau_\epsilon) = \dfrac{\lambda}{2}\tau_\epsilon^{-3/2}e^{-\lambda \tau_\epsilon^{-1/2}},\qquad\text{where }\lambda = -\dfrac{\log(0.5)}{3}\text{ is such that }\pi(\sigma_\epsilon>3) = 0.5 \]
Notice that \(\sigma_u\) and \(\sigma_\epsilon\) are \(\text{Exp}(\lambda)\).
hyper.ar1 = list(theta1 = list(prior = "pc.prec", param = c(0.02, 0.5)), #sigma_u
theta2 = list(prior = "pc.cor1", param = c(0.9, 0.5))) #rho
hyper.family = list(theta = list(prior = "pc.prec", param = c(3, 0.5))) #sigma_epsilon
formula2 <- y ~ f(t,model = 'ar1', hyper = hyper.ar1)
model3 <- inla(formula = formula2,
data = df,
family="gaussian",
control.predictor = list(compute=TRUE),
control.family = list(hyper = hyper.family))
summary(model3)## Time used:
## Pre = 0.149, Running = 0.194, Post = 0.0109, Total = 0.354
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## (Intercept) 8.64 0.305 8.034 8.642 9.239 8.642 0
##
## Random effects:
## Name Model
## t AR1 model
##
## Model hyperparameters:
## mean sd 0.025quant 0.5quant
## Precision for the Gaussian observations 2.20e+04 2.42e+04 1452.689 1.44e+04
## Precision for t 6.12e-01 7.70e-02 0.472 6.09e-01
## Rho for t 7.84e-01 2.80e-02 0.726 7.86e-01
## 0.975quant mode
## Precision for the Gaussian observations 8.64e+04 3952.496
## Precision for t 7.76e-01 0.603
## Rho for t 8.36e-01 0.788
##
## Marginal log-Likelihood: -236.96
## is computed
## Posterior summaries for the linear predictor and the fitted values are computed
## (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')